2南方医科大学南方医院病理科, 广州, 510515
作者 通讯作者
计算分子生物学, 2013 年, 第 2 卷, 第 5 篇
收稿日期: 2013年05月21日 接受日期: 2013年05月21日 发表日期: 2013年05月21日
引用格式(中文):
齐鲁, 丁彦青, 2013, 大肠癌转移相关基因表达调控的生物信息学分析, 基因组学与应用生物学, 32(1): 83-90
引用格式(英文):
Qi L., and Ding Y.Q., 2013, Bioinformatics Analysis of Gene Expression Regulation Associated with Colorectal Cancer Metastasis, Jiyinzuxue Yu Yingyong Shengwuxue (Genomics and Applied Biology), 32(1): 83-90
大肠癌是常见的消化道恶性肿瘤,在中国呈逐年上升的趋势。对大肠癌发生发展转移的研究能够指导临床治疗,对研发新药也有着重要意义。本文通过对表达谱数据进行分析,通过表达谱差异数据进行功 能富集,研究了大肠癌转移前的早期原发肿瘤的转录调控特点以及远隔器官转移后的大肠癌的转录调控特 点,筛选出了部分能够受到多重调控并在转移后肿瘤组织中高表达的关键基因,通过对这些基因相互作用关 系研究,构建出转移后关键基因相互作用调控网络,为大肠癌的治疗提供更多潜在靶点。
蛋白质是基因功能的执行者,基因的转录与翻译对蛋白质的形成有着至关重要的作用,而转录与翻译又受到多种因素的调控和影响。其中转录因子与 DNA 相结合能够调控基因的表达情况,而 miRNA 能够通过结合 RNA 的 3' UTR 区域降解或者抑制 RNA 阻止 RNA 翻译为蛋白质。因此转录水平的调控和翻译水平的调控均能够影响蛋白质最终的形成及产量的多少。
大肠癌是大肠粘膜上皮来源的恶性肿瘤,是消化道恶性肿瘤中的常见类型。因其恶性程度高,易经淋巴道、血道转移,肿瘤晚期常发生远隔器官转移,危及患者生命。转移是大肠癌患者的主要死因。因此研究与大肠癌转移相关重要基因的调控机制,能够为我们提供更多新的药物治疗靶点,以期更好的研究大肠癌的治疗方法。
现有的相关研究主要集中在与大肠癌转移相关的单个基因或 miRNA 的研究上面,或只从差异基因的上下调程度筛选大肠癌转移相关基因,因此忽略了大肠癌转移中各个差异基因的受调控情况。基因受到的转录因子及 miRNA 作用越多,其所受调控越精细, 越能说明此基因的重要性。而在细胞信号的传递过程中,蛋白需与其它蛋白相互作用才能发挥其效能,因此多个重要蛋白构成的相互作用网络,成为细胞信号转导的关键部分。
本文主要通过对大肠癌转移相关差异表达基因的调控情况进行分析和筛选,进而从差异基因中筛选出受调控最多的基因并构建由这些差异基因所构成的调控网络,这些基因在大肠癌转移相关的细胞信号传递中起着关键作用,因此能够成为抑制或者治疗大肠癌转移的重要治疗靶点。
1结果与分析
1.1 miRNA 对大肠癌转移相关差异表达基因的调控
大肠癌 T1、T2 期的基因表达情况与 M1 期有明显的差异,为了更好地挖掘出差异表达基因所包含的信息,通过 GSEA (Damian and Gorfine, 2004)工具对差异表达基因进行不同功能集的基因富集。用 GSEA 通过 miRNA 基因集对此表达谱数据中具有相同 miRNA 靶点的高表达基因进行富集,基因集版本为 V3.0。数据库中的 221 个 miRNA 基因集中有 202 个 可用,并且有 192 个基因集在 M1 表型的表达谱数据中表达上调,只有 10 个基因集在 T1T2 表型的表达谱数据中上调,且假阳性率较高,因此,本文只研究在 M1 表型中上调的基因集。M1 表型的基因集结果通过标准化富集分数进行排序,并且权衡了 P 值及错误发现率(FDR),筛选得到几个在转移后 M1 期的大肠癌高表达相关基因中具有的共同 miRNA 作用靶序列。这些靶序列分别为 MIR-369-3P、MIR-105、MIR-484、 MIR-518B、MIR-518C、MIR-518D、MIR-18A、MIR- 18B、MIR-200B、MIR-200C、MIR-429、MIR-520G、 MIR-520H、MIR-126、MIR-143、MIR-29A、MIR-29B、 MIR-29C、MIR-19A、MIR-19B 能够结合的 3' UTR 的靶序列。每种 miRNA 的靶序列都对应着许多与 M1 期转移相关高表达基因,每种高表达基因的 3' UTR 区域都含有相应上述 miRNA 作用靶序列。虽然每个 miRNA 作用位点集内包含着许多基因,但只有部分在 M1 中的高表达的基因对此功能集起富 集作用(图1)。由于同样的调控序列可能对应着不同的 miRNA,因此上述筛选出的 miRNA 作用位点有11类。将上述 11 类功能集当中起富集作用的基因挑选出来,计算出现频率,研究共调控机制。结果 KLF12 基因的 3' UTR 区域存在上述 11 类中的 6 类 miRNA 作用位点,说明这六类 miRNA 能够同时调控高表达 的 KLF12 基因。同理,具有五类 miRNA 作用位点的基因为 PCDHA9、PCDHA3、PCDHA10、PCDHA5、QKI 、PRICKLE2、ATXN1。具有四类的为 RFX4、OGT、CRE B5、IGF1、TRIM33。对上述 13 种蛋白运用 Visant (Hu et al., 2008)工具构建蛋白质相互作用网络图(图2),可以看出,AXIN1 可以将其中的 7 个蛋白连接起来,AXIN1 为 Wnt 信号通路中的轴蛋白(Biechele et al., 2012),因此,AXIN1 在大肠癌转移相关基因的调控上起着重要作用。
图1 GSEA富集图 |
图2 13个蛋白与蛋白间相互作用网络 |
为了更好的分析筛选出的 13 个关键蛋白在大肠癌转移后分子网络中的关系,将上述 13 种蛋白的相互作用蛋白每个节点进一步展开,得到一个更复杂的蛋白相互作用网络图,然后计算能把上述蛋白都关联 起来的中间节点的蛋白数量共 364 个。上述 13 种蛋白可能在大肠癌转移起着重要作用,但每种蛋白并不是单独能够完成相应的功能,而是通过与它相互作用 的蛋白构成多级网络,将生物信号进行级联传递的。 所以 364 个蛋白作为桥梁,把这 13 种蛋白联系起来, 构成一个复杂的网络,而这网络,可能为大肠癌转移 发生发展相关的信号传递网络的一部分。因此对这些蛋白进行 Gene Set Analysis Toolkit V2 工 具 中 的 Wikipathways (Kelder et al., 2012)通路分析(表1)可以发现,雄激素受体信号通路,TGF-β 信号通路,EGFR1 信号通路等均有许多中间节点蛋白参与,因此,雄激 素受体信号等通路与大肠癌的转移可能有着密切联系,因此也可以认为激素水平的变化对大肠癌的转移可能有着重要影响。对这 364 个蛋白共同转录因子通过 Gene Set Analysis Toolkit V2 工具进行富集分析,结果 364 个基因当中有 102 个基因上游调控区域有 SP1 转录因子结合位点, 有 86 个基因有 MAZ 结合位点,有 61 个基因有 PAX4 结合位点,有 69 个基因有 LEF1 结合位点,有 72 个基因有 E12 结合位点,有 65 个基因有 FOXO4 结合位点,有 57 个 基因有 AP4 结合位点,有 57 个基因有 NFAT 结合位点,有 39 个基因有 MYOD 结合位点,这些转录因子结合位点调控着这 364 个基因构成的网络,可能与大肠癌转移密切相关。
表1维基通路富集分析 |
1.2转录因子对大肠癌转移相关差异表达基因的调控
用 GSEA 工具通过转录因子结合位点基因集对 前面已归为 T1T2,M1 两类的大肠癌表达谱数据进行共同转录因子结合位点富集。在 615 个转录因子基因集中,有 576 个可用,其中有 486 个基因集在 M1 表 型中表达上调,只有 90 个基因集在 T1T2 表型中表 达上调,可见转移后的肿瘤高表达基因需要较多的转录因子进行精密调控。在 M1 表型表达上调的 486 个 转录因子结合位点基因集中,通过标准化富集分数, 权衡 P 值及错误发现率,经筛选得到前几个转录因子结合位点分别是:EVI1、OCT1、FOXO4、SMAD3、AM L、 PAX2、TAL1ALPHAE47、FREAC2、PAX4、AM-EF2、 RREB1、CEBPA、TEF、SRY、GATA1、FXR 等的转录因子的结合位点,说明在转移后的大肠癌组织中,表达升高的基因均含有上述转录因子结合位点,这些对应的转录因子对大肠癌转移起着重要作用。而部分基因的转录调控区域同时含有这些转录因子结合位点,构成了组合调控,调控机制越复杂,转录因子越多,说明这些基因表达调控越精细,这部分基因在转移后的大肠癌组织中发挥的功能越重要。在此 16 个转录因子结合位点基因集中,提取对相应基因集起着富集作用的基因,计算得出同时含有 4 个以上此转录因子结合位点的基因共有 41 个,其中有 HOXC4,PURA 基因同时含有 9 个此转录因子结合位点。因此,这些基因在转移后相关转录因子中,受到多个相关转录因子组合调控,能够对大肠癌转移的发生发展起到关键作用。对这 41 个有着大于 4 个相关转录因子结合位点的基因进行 GO 分析,得出其中有 22 个基因能够与核酸结合,其中有 19 个具有转录活性,为转录因子(图3),说明这些转录因子同时也调控着其他基因。对这 41 个基因对应的蛋白构建相互作用网络,将每个与这 41 个蛋白相连的蛋白进一步展开,然后计算与这 41 个蛋白相连的中间节点蛋白,共有 1 100 个,也就是说有 1 100 个蛋白把这 41 个蛋白连接起来构成复杂的相互作用网络,对这 1 100 个蛋白进行 Gene Set Analysis Toolkit V2 工具中的 Wikipathways 通路分析(表2)。可以发现,TGF-β 通路,雄激素受体信号 通路等与这个网络相关,且与之前 miRNA 调控基因所构成的网络一致,因此进一步说明了大肠癌的转移与这几个信号通路相关。雄激素为固醇类激素,可直接通过细胞膜与细胞核内转录因子结合,调控转录。 因此,雄激素调控转录因子,转录因子又通过相互作用网络将信号进一步传递,调控基因的表达。对这 1 100 个蛋白通过 Gene Set Analysis Toolkit V2 进行共同转录因子进行分析,结果有 277 个基因有 SP1 结合位点,有 240 个基因有 MAZ 结合位点,有 185 个基因有 LEF1 结合位点,有 175 个基因有 FOXO4 结合位点,有 167 个基因有 NFAT 结合位点,有 186 个基因有 E12 结合位点,有 138 个基因有 PAX4 结合位 点,有 119 个基因有 NFY 结合位点,有 101 个基因有 FRAC2 结合位点。可以发现,大部分转录因子与前面 miRNA 所调控基因构成网络的共同转录因子重叠。 说明这些转录因子在大肠癌转移相关信号调控网络中起着重要作用。
图3分子功能与细胞组分条形图 |
表2维基通路富集分析 |
综合前面转移后大肠癌中 miRNA 所调控的基因 以及转录因子所调控的基因,可以发现基因 KLF12 以 及 CREB5 不仅在转录调控区域有很多相关转录因子结合位点,而且在转录后的 3' UTR 区域有着很多相关 miRNA 结合位点,因此认为这两个基因能够同时被这 些转录因子以及 miRNA 所调控(图4)。因为 KLF12 与 CREB5 本身为转录因子,因此也可以调控其他基因转录,形成级联反应,构建 KLF12 与 CREB5 相互作用调 控网络图并将与他们相互作用蛋白进行展开,计算能够将 KLF12 与 CREB5 关联起来的若干个蛋白,并构成新的网络如(图5)。可以看出在 KLF12 与 CREB5 之间有 28 个与这两个转录因子直接或间接相互作用的蛋白,这些蛋白将其联系在一起。在此对这 28 个蛋白 进行 GO 分析,28 个蛋白里面具有转录调控活性的蛋白有 26 个,因此可以说明,这两个蛋白是通过这 26 个 转录因子连接起来的,而转录因子相互结合成复合体可调控其他基因表达,因此这些基因与大肠癌转移密切相关。
图4 KLF12 和 CREB5 的 miRNA 和转录因子相互作用网络 |
图5 KLF12 和 CREB5 蛋白与蛋白相互作用网络 |
在前面的 GSEA 分析中, 有 90 个基因集在 T1T2 表型中表达上调。说明部分基因在转移前的原发肿瘤早期起着重要作用,通过标准化富集分数,权衡 P 值及错误发现率。筛选得到最相关的几个转录因子结合位点基因集为 E2F_Q6、E2F1_Q3、E2F1_Q6_01、E2F_ Q4、E2F1_Q6、E2F_03、E2F1DP2_01、E2F_02、E2F1DP1_01、E2F1DP2_01、E2F4DP2_01、E2F4DP1_01。可以发现,全部转录因子结合位点均为 E2F 家族,或者 E2F 与 DP 家族的二聚体,其中主要为 E2F1。可见在转移前的原发肿瘤早期,E2F 家族转录因子对肿瘤的发生发展起着极为关键的作用,提取对相应基因集起着富集作用的基因,经计算有 18 个基因均含有上述12 种转录因子结合位点,有 62 个基因同时含有 8 个 以上的上述转录因子结合位点,说明这些基因在原发肿瘤早期起着关键作用。因此将这 62 个基因进行 Gene Set Analysis Toolkit V2 工具中的 Wikipathways 通路分析(表3),结果可以看出,在转移前原发肿瘤早期,E2F 家族主要转录调控的高表达基因主要参与了细胞代谢例如 DNA 复制,核苷酸合成,染色体重组等生物进程,因此,说明在原发肿瘤早期主要以细胞增殖,代谢增强为主要表现。
表3维基通路富集分析 |
2讨论
大肠癌的发生发展转移是由多因素多步骤共同作用的结果,发生远隔器官转移的大肠癌与未发生转移的原位大肠癌基因表达有着明显差异,对这些差异表达基因的研究,能够使我们更好的了解大肠癌转移的分子机制。本研究不仅限于对大肠癌 T1T2 期与 M1 期差异度高的基因进行分析,而是通过对差异度高的基因的调控区域的 miRNA 作用位点和转录因子结合位点进行功能富集,筛选出差异基因共有的 miRNA 作用位点和转录因子结合位点,并将含有这些相关的作用位点较多的差异基因进行筛选,因此筛选出的基因受到众多相关的 miRNA 和转录因子的调控,功能越重要的基因受到的调控越精细,因此这些基因对大肠癌的转移有着重要作用。对这些基因构建相互作用网络,再将网络中的节点蛋白共有的转录因子结合位点进行富集,可以筛选出对此蛋白网络起着重要作用的转录因子。因此,本研究筛选出 2 个含有这些作用位点最多的差异表达基因 KLF12 和 CRE B5。这两个基因本身为转录因子,连接他们的相互作用网络中大部分也为转录因子,因此这些转录因子群相互连接,组合调控其他基因的表达,显示出其对大肠癌转移相关基因表达调控的影响力。KLF12 基因已被证实在胃癌中表达升高,并且高表达 KLF12 能够增强胃癌的侵袭和转移,证明其在低分化胃癌中能 够促进肿瘤的发生发展(Nakamura et al., 2009)。CREB5 为 CRE 依赖的转录因子,它是唯一一个能够连接 cAMP 信号通路和 TPA 信号通路的关键蛋白(Zu et al., 1993),其与大肠癌发生发展转移的关系有待于进一 步研究。本文构建了蛋白质相互作用网络,并且发现 调控网络中蛋白质所对应基因的上游调控区域转录因子结合位点较为一致,因此筛选出对此蛋白网络起着重要作用的转录因子包括 SP1、MAZ、LEF1、FOXO4、 NFAT、E12、PAX4 等。其中 SP1 能够调控多种肿瘤相关基因的表达,参与了肿瘤的增值,分化和凋亡,并能够调控雄激素受体以及 TGF-β (Sankpal et al., 20 11),这与前面通路富集的结果一致。对关键基因相互作用网络的中间蛋白进行通路富集分析,可以看出 TGF-β 通路,雄激素受体信号通路等与网络中的蛋白密切相关,说明与大肠癌转移相关的蛋白相互作用网络中这些信号通路参与了重要作用,TGF-β 通路与肿瘤的 EMT相关,而雄激素受体信号通路与体内激素水平密切相关。E2F1 转录因子结合位点在 T1T2 期的大肠癌高表达基因转录调控区域中都存在并富集,说明 E2F1 基因在早期大肠癌对肿瘤的生长起着重要作用。E2F1 在细胞周期中参与了 DNA 的合成、 修复,细胞的 G1 期到 S 期的调定点,以及有丝分裂和凋亡等功能(Zacharatos et al., 2004),这与前面通路富集结果一致。 高 表 达 E2F1 能够促进肿瘤的生长 (Zacharatos et al., 2004),因此,在大肠癌转移前, E2F1 能够对肿瘤的增殖起重要作用。大肠癌的发生发展转移是一个由多组基因,多条信号通路共同作用,相互调控的结果,其涉及的基因众多,调控关系复杂,在此通过生物信息学的方法,能够很好的找到表达调控的关键节点,对指导实验研究有着重要意义。本文 寻找到了多个基因表达调控的关键节点,对这些关键节点上进行干预,或许能够减缓甚至抑制大肠癌的发生发展转移,为临床的有效治疗提供更好的治疗靶点。
3材料与方法
3.1数据来源
本文选取的是来自 NCBI 网站的 GEO 数据库中芯片数据 GSE2109 系列的大肠癌数据进行分析,此系列本身包含 2 158 个各类肿瘤的数据并包含各种临床信息。由于本研究主要集中于大肠癌,因此在这 2 158例表达谱数据中挑选大肠癌数据并同时包含 TNM 分期数据 的表达谱数据进行分析,可挑选出 343 例大肠癌数据,为了更好的分析未发生转移的肿瘤与已发生远隔器官转移肿瘤的基因表达差异情况,从 343 例大肠癌表达谱数据中挑选出归类于 T1、T2 和 M1 分期的大肠癌表达谱数据共 121 例进行分析,将 TNM 分期中属于 T1、T2 的归为一类,M1 的归为一类。
3.2分析方法
本文主要的实验方法是通过 GSEA 筛选出与大肠癌转移相关的 miRNA 和转录因子功能集,由于每个功能集包含很多差异表达基因,因此对筛选出的不同功能集中的多个差异表达基因取交集,取交集后的差异表达基因通过 Gene Set Analysis Toolkit V2 工具进行 GO 分析,归类,分析其功能类。然后对这些差异表达基因所对应的蛋白通过 visant 工具构建蛋白质相互作用网络,提取出除了差异基因对应蛋白以外的相互作用蛋白,再对这些作用蛋白通过 Gene Set Analysis Toolkit V2 工具进行 Wikipathways 通路富集和转录因子富集。
作者贡献
第一作者齐鲁对本论文进行设计、数据搜集、分析以及论文撰写;通讯作者丁彦青教授对本论文进行指导及审阅。
致谢
感谢评审专家对论文的评审和修改建议。
Biechele T.L., Kulikauskas R.M., Toroni R.A., Lucero O.M., Swift R. D., James R.G., Robin N.C., Dawson D.W., Moon R.T., and Chien A.J., 2012, Wnt/beta-catenin signaling and AXIN1 regu- late apoptosis triggered by inhibition of the mutant kinase BRAFV600E in human melanoma, Sci. Signal., 5(206): ra3
Damian D., and Gorfine M., 2004, Statistical concerns about the GSEA procedure, Nat. Genet., 36(7): 663
Hu Z., Snitkin E.S., and DeLisi C., 2008, VisANT: An integra- tive framework for networks in systems biology, Brief Bioinform, 9(4): 317-325
Kelder T., van Iersel M.P., Hanspers K., Kutmon1 M., Conklin B.R., Chris T., Evelo C.T., and Pico A.R., 2012, WikiPath- ways: Building research communities on biological path- ways, Nucleic Acids Res., 40(D1): 1301-1307
Nakamura Y., Migita T., Hosoda F., Okada N., Gotoh M., Arai Y., Fukushima M., Ohki M., Miyata S., Takeuchi K., Imoto I., Katai H., Yamaguchi T., Inazawa J., Hirohashi S., Ishikawa Y., and Shibata T., 2009, Krüppel-like factor 12 plays a significant role in poorly differentiated gastric cancer progression, Int. J. Cancer, 125(8): 1859-1867
Sankpal U.T., Goodison S., Abdelrahim M., and Basha R., 2011, Targeting sp1 transcription factors in prostate cancer thera- py, Med. Chem., 7(5): 518-525
Zacharatos P., Kotsinas A., Evangelou K., Karakaidos P., Vassil- iou L.V., Rezaei N., Kyroudi A., Kittas C., Patsouris E., Pa-pavassiliou A.G., and Gorgoulis V.G., 2004, Distinct ex- pression patterns of the transcription factor E2F-1 in rela- tion to tumour growth parameters in common human carci-nomas, J. Pathol., 203(3): 744-753
Zu Y.L., Maekawa T., Nomura N., Nakata T., and Ishii S., 1993, Reg- ulation of trans-activating capacity of CRE-BPa by phorbol es- ter tumor promoter TPA, Oncogene, 8(10): 2749-2758